https://femiguez.github.io/nlraa-docs/nlraa-Oddi-LFMC.html
https://onlinelibrary.wiley.com/doi/pdfdirect/10.1002/ece3.5543
#install.packages("emmeans")
library(nlme)
library(nlraa)
library(ggplot2)
library(emmeans) # emmeans::contrast(), emmeans::emmeans()
data(lfmc)
lfmc.gd <- groupedData(lfmc ~ time | group, data = lfmc, order.groups = FALSE)
lfmc.gd$plot = with(lfmc.gd, factor(plot, levels =c("4", "5", "6", "1", "2", "3")))
head(lfmc.gd)
Grouped Data: lfmc ~ time | group
leaf.type time plot site lfmc group
1 M. spinosum 1 1 P 238.4 SM plot 1
2 M. spinosum 1 1 P 269.8 SM plot 1
3 M. spinosum 1 1 P 325.2 SM plot 1
4 S. bracteolactus 1 1 P 326.7 SS plot 1
5 S. bracteolactus 1 1 P 257.3 SS plot 1
6 S. bracteolactus 1 1 P 279.3 SS plot 1
ggplot(lfmc, aes(time, lfmc)) +
geom_point() +
geom_smooth(method = "loess") +
facet_wrap(~leaf.type) +
xlab("Time (days)") +
ylab("LMFC (%)")
`geom_smooth()` using formula 'y ~ x'
nlsL <- nlsList(lfmc ~ SSdlf(time, A, w, m, s), data = lfmc.gd)
Warning: 3 errors caught in nls(model, data = data, control = controlvals). The error messages and their frequencies are
step factor 0.000488281 reduced below 'minFactor' of 0.000976562
1
number of iterations exceeded maximum of 50
2
coef(nlsL)
plot(intervals(nlsL))
nl.re <- nlme(nlsL, control = nlmeControl(maxIter = 5000, msMaxIter = 1500))
Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: function evaluation limit reached without convergence (9)
nl.re.1 <- update(nl.re, random = pdDiag(A + w + m + s ~ 1))
fxf <- fixef(nl.re.1)
fxf
A w m s
193.45254 21.45745 31.48273 -22.59708
# nl.me <- update(nl.re.1, fixed = list(A + w + m + s ~ leaf.type),
# start = c(A = fxf[1],0,0,0,
# w = fxf[2],0,0,0,
# m = fxf[3],0,0,0,
# s = fxf[4],0,0,0))
nl.re.2 <- update(nl.re.1, random = pdDiag(A + w + m ~ 1))
fxf <- fixef(nl.re.2)
fxf
A w m s
193.46654 21.43199 31.48852 -22.60602
nl.me.1 <- update(
nl.re.2,
fixed = list(A + w + m + s ~ leaf.type),
start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3],0,0,0, s=fxf[4],0,0,0))
nl.me.1
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -1034.375
Fixed: list(A + w + m + s ~ leaf.type)
A.(Intercept) A.leaf.typeGrass E
23.135414 53.294908
A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
296.620596 263.313118
w.(Intercept) w.leaf.typeGrass E
50.291577 -34.313924
w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
-26.732287 2.135665
m.(Intercept) m.leaf.typeGrass E
51.533983 -21.920404
m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
-20.647403 -18.170359
s.(Intercept) s.leaf.typeGrass E
17.605794 -25.926031
s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus
-43.529450 -33.666226
Random effects:
Formula: list(A ~ 1, w ~ 1, m ~ 1)
Level: group
Structure: Diagonal
A.(Intercept) w.(Intercept) m.(Intercept) Residual
StdDev: 8.112813 0.001245815 2.468334 15.38607
Number of Observations: 247
Number of Groups: 12
nl.re.2
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -1070.648
Fixed: list(A ~ 1, w ~ 1, m ~ 1, s ~ 1)
A w m s
193.46654 21.43199 31.48852 -22.60602
Random effects:
Formula: list(A ~ 1, w ~ 1, m ~ 1)
Level: group
Structure: Diagonal
A w m Residual
StdDev: 119.3792 10.94165 8.927279 15.26574
Number of Observations: 247
Number of Groups: 12
AIC(nl.re.2, nl.me.1)
IC_tab(nl.re.2, nl.me.1)
nl.me.2 <- update(nl.me.1, random = pdDiag(A + m ~ 1))
IC_tab(nl.me.1, nl.me.2)
nl.me.3 <- update(nl.me.2, random = pdDiag(A ~ 1))
IC_tab(nl.me.2, nl.me.3)
plot(nl.me.3)
rse <- residuals(nl.me.3) / sigma(nl.me.3)
hist(residuals(nl.me.3, type = "normalized"))
#par(pty = "s")
qqnorm(rse)
qqline(rse)
nl.me.3.vm <- update(nl.me.3, weights = varIdent(form = ~1 | leaf.type))
IC_tab(nl.me.3, nl.me.3.vm)
plot(nl.me.3.vm)
hist(residuals(nl.me.3.vm, type = "normalized"), main = "", xlab = "Residuals")
qqnorm(resid(nl.me.3.vm, type = "normalized"))
qqline(resid(nl.me.3.vm, type = "normalized"))
plot(ACF(nl.me.3.vm), alpha = 0.01)
nl.me.ml <- update(nl.me.3.vm, method = "ML")
nl.me.ml
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -965.5628
Fixed: list(A + w + m + s ~ leaf.type)
A.(Intercept) A.leaf.typeGrass E
50.1522850 28.9553860
A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
346.5816365 231.1569668
w.(Intercept) w.leaf.typeGrass E
24.6693608 -9.4837876
w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
-35.8327959 29.7533014
m.(Intercept) m.leaf.typeGrass E
48.9613646 -19.9780200
m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
-26.2920797 -14.9909275
s.(Intercept) s.leaf.typeGrass E
-16.2294126 6.3520529
s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus
-21.0664461 0.9389322
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.412262 6.824811
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8547836 3.1456077 3.0262123
Number of Observations: 247
Number of Groups: 12
nl.me.ml.1 <- update(nl.me.ml,
fixed = list(A + w + m ~ leaf.type, s ~ 1),
start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], 0, 0, 0, s = fxf[4])
)
nl.me.ml.1
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -969.2617
Fixed: list(A + w + m ~ leaf.type, s ~ 1)
A.(Intercept) A.leaf.typeGrass E
49.94562 39.60979
A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
219.45652 232.11909
w.(Intercept) w.leaf.typeGrass E
25.24501 -14.34293
w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
35.36121 28.70388
m.(Intercept) m.leaf.typeGrass E
48.30390 -21.59477
m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
-15.57765 -14.37454
s
-15.44291
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.612472 6.825604
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8714808 3.2750695 3.0238813
Number of Observations: 247
Number of Groups: 12
IC_tab(nl.me.ml, nl.me.ml.1)
nl.me.ml.2 <- update(nl.me.ml, fixed = list(A + w ~ leaf.type, m + s ~ 1),
start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3], s=fxf[4]))
nl.me.ml.2
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -972.7465
Fixed: list(A + w ~ leaf.type, m + s ~ 1)
A.(Intercept) A.leaf.typeGrass E
54.25730 30.93396
A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
223.27914 240.31543
w.(Intercept) w.leaf.typeGrass E
29.05803 -20.68663
w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
31.68646 26.99003
m s
30.91655 -16.15435
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.216113 7.028345
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8820413 3.1693967 2.9643305
Number of Observations: 247
Number of Groups: 12
nl.me.ml.3 <- update(nl.me.ml, fixed = list(A ~ leaf.type, w + m + s ~ 1),
start = c(A=fxf[1],0,0,0, w=fxf[2], m=fxf[3], s=fxf[4]))
nl.me.ml.3
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -1017.91
Fixed: list(A ~ leaf.type, w + m + s ~ 1)
A.(Intercept) A.leaf.typeGrass E
63.92214 14.27256
A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
251.99464 265.42794
w m
17.33007 32.65667
s
-25.16304
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 8.13893 7.795662
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.000000 1.584893 2.813972 2.731419
Number of Observations: 247
Number of Groups: 12
nl.me.ml.4 <- update(nl.me.ml, fixed = list(w ~ leaf.type, A + m + s ~ 1),
start = c(A=fxf[1], w=fxf[2],0,0,0, m=fxf[3], s=fxf[4]))
nl.me.ml.4
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -999.8224
Fixed: list(w ~ leaf.type, A + m + s ~ 1)
w.(Intercept) w.leaf.typeGrass E
29.16028 -20.38030
w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
33.19422 28.72473
A m
175.78134 31.08599
s
-15.70300
Random effects:
Formula: A ~ 1 | group
A Residual
StdDev: 108.0655 7.042353
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8770644 3.1888994 2.9304282
Number of Observations: 247
Number of Groups: 12
IC_tab(nl.me.ml.1, nl.me.ml.2, nl.me.ml.3, nl.me.ml.4)
M1 <- update(nl.me.ml.2, method = "REML")
summary(M1)
Nonlinear mixed-effects model fit by REML
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.408521 7.162899
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8851762 3.1796110 2.9647390
Fixed effects: list(A + w ~ leaf.type, m + s ~ 1)
Correlation:
A.(In) A.l.GE A..M.s A..S.b w.(In) w.l.GE w..M.s
A.leaf.typeGrass E -0.527
A.leaf.typeM. spinosum -0.146 0.535
A.leaf.typeS. bracteolactus -0.115 0.537 0.746
w.(Intercept) -0.217 -0.034 -0.200 -0.216
w.leaf.typeGrass E -0.035 -0.283 -0.382 -0.399 -0.258
w.leaf.typeM. spinosum -0.118 -0.227 -0.547 -0.454 0.157 0.573
w.leaf.typeS. bracteolactus -0.129 -0.241 -0.462 -0.575 0.188 0.601 0.640
m -0.217 -0.324 -0.633 -0.666 0.092 0.145 0.161
s -0.246 -0.354 -0.710 -0.748 0.416 0.575 0.702
w..S.b m
A.leaf.typeGrass E
A.leaf.typeM. spinosum
A.leaf.typeS. bracteolactus
w.(Intercept)
w.leaf.typeGrass E
w.leaf.typeM. spinosum
w.leaf.typeS. bracteolactus
m 0.172
s 0.752 0.544
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.27999716 -0.59657157 -0.02042043 0.57273326 3.12017921
Number of Observations: 247
Number of Groups: 12
fixef(M1)
A.(Intercept) A.leaf.typeGrass E
54.25350 30.92293
A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
223.24504 240.27654
w.(Intercept) w.leaf.typeGrass E
29.05581 -20.68938
w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
31.67297 26.97508
m s
30.92881 -16.15406
ranef(M1)
A.(Intercept)
GW plot 4 -1.915443
GW plot 5 -1.923198
GW plot 6 3.838641
GE plot 1 15.652459
SM plot 1 6.967303
SS plot 1 9.166925
GE plot 2 -11.066197
SM plot 2 -5.358910
SS plot 2 2.442781
GE plot 3 -4.586263
SM plot 3 -1.608393
SS plot 3 -11.609707
intervals(M1)
Approximate 95% confidence intervals
Fixed effects:
lower est. upper
A.(Intercept) 42.45370 54.25350 66.05331
A.leaf.typeGrass E 13.51301 30.92293 48.33286
A.leaf.typeM. spinosum 191.88318 223.24504 254.60691
A.leaf.typeS. bracteolactus 207.05884 240.27654 273.49423
w.(Intercept) 25.68039 29.05581 32.43122
w.leaf.typeGrass E -25.79788 -20.68938 -15.58088
w.leaf.typeM. spinosum 16.39521 31.67297 46.95073
w.leaf.typeS. bracteolactus 11.07083 26.97508 42.87934
m 26.85931 30.92881 34.99831
s -20.81510 -16.15406 -11.49302
Random Effects:
Level: group
lower est. upper
sd(A.(Intercept)) 5.035417 9.408521 17.57953
Variance function:
lower est. upper
Grass E 0.6772617 0.8851762 1.156919
M. spinosum 2.4528838 3.1796110 4.121649
S. bracteolactus 2.2864085 2.9647390 3.844316
Within-group standard error:
lower est. upper
5.970704 7.162899 8.593144
# 2.2.1 - Overall predictions (Table 3) ----
Agw <- fixef(M1)[1] # "A" for grasses in the W site (GW)
Age <- fixef(M1)[1]+fixef(M1)[2] # "A" for grasses in the E site (GE)
Asm <- fixef(M1)[1]+fixef(M1)[3] # "A" for Mullinum spinosum (SM)
Ass <- fixef(M1)[1]+fixef(M1)[4] # "A" for Senecio filaginoides (SS)
Wgw <- fixef(M1)[5] # "w" for grasses in the W site (GW)
Wge <- fixef(M1)[5]+fixef(M1)[6] # "w" for grasses in the E site (GE)
Wsm <- fixef(M1)[5]+fixef(M1)[7] # "w" for Mullinum spinosum (SM)
Wss <- fixef(M1)[5]+fixef(M1)[8] # "w" for Senecio filaginoides (SS)
M <- fixef(M1)[9] # "m" for all the leaf types
S <- fixef(M1)[10] # "s" for all the leaf types
GW <- subset(lfmc, leaf.type == "Grass W")
GE <- subset(lfmc, leaf.type == "Grass E")
SM <- subset(lfmc, leaf.type == "M. spinosum")
SS <- subset(lfmc, leaf.type == "S. bracteolactus")
par(mfrow=c(2,2), cex=0.75)
# Grasses W
plot(GW$time, GW$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses W site", col="darkgreen")
curve(((Agw-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=2, lty=1, add=T, col="darkgreen")
# Grasses E
plot(GE$time, GE$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses E site", col="green")
curve(((Age-Wge)/(1 + exp((M-x)/S))+Wge), lwd=2, lty=1, add=T, col="green")
# M. Spinosum
plot(SM$time, SM$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="M. spinosum", font.main=4, col="darkorange")
curve(((Asm-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=2, lty=1, add=T, col="darkorange")
# S. filaginoides
plot(SS$time, SS$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="S. bracteolactus", font.main=4, col="darkred")
curve(((Ass-Wss)/(1 + exp((M-x)/S))+Wss), lwd=2, lty=1, add=T, col="darkred")
# 2.2.2 - Predictions at plot level (Table 3) ----
Agw.p4 <- fixef(M1)[1]+ranef(M1)[1,1] # "A" for grasses in the plot 4 (GW)
Agw.p5 <- fixef(M1)[1]+ranef(M1)[2,1] # "A" for grasses in the plot 5 (GW)
Agw.p6 <- fixef(M1)[1]+ranef(M1)[3,1] # "A" for grasses in the plot 6 (GW)
Age.p1 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[4,1] # "A" for grasses in the plot 1 (GW)
Age.p2 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[7,1] # "A" for grasses in the plot 2 (GW)
Age.p3 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[10,1] # "A" for grasses in the plot 3 (GW)
Asm.p1 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[5,1] # "A" for M. spinosum in the plot 1 (SM)
Asm.p2 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[8,1] # "A" for M. spinosum in the plot 2 (SM)
Asm.p3 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[11,1] # "A" for M. spinosum in the plot 3 (SM)
Ass.p1 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[6,1] # "A" for S. filaginoides in the plot 1 (SM)
Ass.p2 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[9,1] # "A" for S. filaginoides in the plot 2 (SM)
Ass.p3 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[12,1] # "A" for S. filaginoides in the plot 3 (SM)
par(mfrow=c(2,2), cex=0.75) # (Figure 4)
# Grasses W site
plot(GW$time, GW$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses W site", col="darkgreen")
curve(((Agw.p4-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 4
curve(((Agw.p5-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 5
curve(((Agw.p6-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 6
# Grasses E site
plot(GE$time, GE$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses E site", col="green")
curve(((Age.p1-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 1
curve(((Age.p2-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 2
curve(((Age.p3-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 3
# M. Spinosum (E site)
plot(SM$time, SM$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="M. spinosum", font.main=4, col="darkorange")
curve(((Asm.p1-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 1
curve(((Asm.p2-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 2
curve(((Asm.p3-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 3
# S. filaginoides (E site)
plot(SS$time, SS$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="S. bracteolactus", font.main=4, col="darkred")
curve(((Ass.p1-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 1
curve(((Ass.p2-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 2
curve(((Ass.p3-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 3
# 2.2.3 - Derivatives (drying speed) ----
lfmc.GW <- expression((Agw-Wgw)/(1 + exp((M-x)/S))+Wgw) # LFMC dynamics for grasses in the W site
ds.GW <- deriv(lfmc.GW, "x", function.arg = TRUE) # Drying speed for grasses in the W site
lfmc.GE <- expression((Age-Wge)/(1 + exp((M-x)/S))+Wge) # LFMC dynamics for grasses in the E site
ds.GE <- deriv(lfmc.GE, "x", function.arg = TRUE) # Drying speed for grasses in the E site
lfmc.SM <- expression((Asm-Wsm)/(1 + exp((M-x)/S))+Wsm) # LFMC dynamics for M. spinosum
ds.SM <- deriv(lfmc.SM, "x", function.arg = TRUE) # Drying speed for M. spinosum
lfmc.SS <- expression((Ass-Wss)/(1 + exp((M-x)/S))+Wss) # LFMC dynamics for S. filaginoides
ds.SS <- deriv(lfmc.SS, "x", function.arg = TRUE) # # Drying speed for S. filaginoides
xvec <- seq(0, 100, length = 1000)
y.ds.GW <- ds.GW(xvec)
y.ds.GE <- ds.GE(xvec)
y.ds.SM <- ds.SM(xvec)
y.ds.SS <- ds.SS(xvec)
par(mfrow=c(2,2), cex=0.75) # (Figure 4)
plot(xvec, y.ds.GW, xlim=c(0,80), ylim=c(0,4),
ylab="Drying speed", xlab="Time (days)", type="n", main="Grasses in the W site")
lines(xvec, -1*(attr(y.ds.GW, "grad")), lwd=2, lty=1, col="darkgreen")
plot(xvec, y.ds.GE, xlim=c(0,80), ylim=c(0,4),
ylab="Drying speed", xlab="Time (days)", type="n", main="Grasses in the E site")
lines(xvec, -1*(attr(y.ds.GE, "grad")), lwd=2, lty=1, col="green")
plot(xvec, y.ds.SM, xlim=c(0,80), ylim=c(0,4),
ylab="Drying speed", xlab="Time (days)", type="n", main="M. spinosum", font.main=4)
lines(xvec, -1*(attr(y.ds.SM, "grad")), lwd=2, lty=1, col="darkorange")
plot(xvec, y.ds.SS, xlim=c(0,80), ylim=c(0,4),
ylab="Drying speed", xlab="Time (days)", type="n", main="S. filaginoides", font.main=4)
lines(xvec, -1*(attr(y.ds.SS, "grad")), lwd=2, lty=1, col="darkred")
# 2.3 ---- Testing differences among leaf types -------------------------------------------
# 2.3.1 - Parameter "A" ----
contrast(emmeans(M1, ~leaf.type, param = "A"), "pairwise")
contrast estimate SE df t.ratio p.value
Grass W - Grass E -30.9 8.84 226 -3.500 0.0031
Grass W - M. spinosum -223.2 15.92 226 -14.027 <.0001
Grass W - S. bracteolactus -240.3 16.86 226 -14.254 <.0001
Grass E - M. spinosum -192.3 13.45 226 -14.294 <.0001
Grass E - S. bracteolactus -209.4 14.22 226 -14.718 <.0001
M. spinosum - S. bracteolactus -17.0 11.71 226 -1.454 0.4671
P value adjustment: tukey method for comparing a family of 4 estimates
# 2.3.2 - Parameter "w" ----
contrast(emmeans(M1, ~leaf.type, param = "w"), "pairwise")
contrast estimate SE df t.ratio p.value
Grass W - Grass E 20.7 2.59 226 7.981 <.0001
Grass W - M. spinosum -31.7 7.75 226 -4.085 0.0004
Grass W - S. bracteolactus -27.0 8.07 226 -3.342 0.0053
Grass E - M. spinosum -52.4 6.62 226 -7.914 <.0001
Grass E - S. bracteolactus -47.7 6.84 226 -6.974 <.0001
M. spinosum - S. bracteolactus 4.7 6.72 226 0.699 0.8974
P value adjustment: tukey method for comparing a family of 4 estimates
# 3.1.1 - Base nonlinear model ----
nl.fe.0 <- nls(lfmc ~ ((A - w) / (1 + exp((m - time)/s))) + w,
start = c(A=15, m=30, s=-17, w=5),
data = lfmc) # here, time is the only predictor variable
b <- coef(nl.fe.0) # coefficients of the base model
# 3.1.2 - Leaf type fixed effect ----
nl.fe.1 <- nls(lfmc ~ ((A[leaf.type] - w[leaf.type])/(1 + exp((m - time)/s))) + w[leaf.type],
start = list(A=rep(b[1], 4), m=25, s=-10, w=rep(b[4], 4)),
data = lfmc) # and fit is made using the base model's coefficients as start values
b1 <- coef(nl.fe.1) # coefficients of the model with leaf type and time as predictor variables
# 3.1.3 - Plot fixed effect ----
nl.fe.2 <- nls(lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group],
start = list(A=rep(c(b1[1],b1[2],b1[3],b1[4]),3), m=25, s=-10, w=rep(c(b1[7],b1[8],b1[9],b1[10]),3)),
data = lfmc) # the model is fit using "b1" as the start values
IC_tab(nl.fe.0, nl.fe.1, nl.fe.2)
# Residual checking:
par(mfrow=c(1,1))
plot(nl.fe.2) # heterocedasticity in the residuals
hist(resid(nl.fe.2), main="", xlab="Residuals") # slightly skew distribution
qqnorm(resid(nl.fe.2))
qqline(resid(nl.fe.2))
M2 <- nl.fe.2
summary(M2)
Formula: lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group]
Parameters:
Estimate Std. Error t value Pr(>|t|)
A1 53.683 8.602 6.240 2.20e-09 ***
A2 54.282 8.636 6.285 1.72e-09 ***
A3 61.906 8.881 6.971 3.59e-11 ***
A4 111.561 13.291 8.394 5.65e-15 ***
A5 314.131 24.840 12.646 < 2e-16 ***
A6 333.831 26.640 12.531 < 2e-16 ***
A7 77.465 10.997 7.044 2.34e-11 ***
A8 298.161 24.917 11.966 < 2e-16 ***
A9 321.433 25.765 12.475 < 2e-16 ***
A10 87.553 11.744 7.455 2.01e-12 ***
A11 279.437 20.829 13.416 < 2e-16 ***
A12 292.443 24.178 12.095 < 2e-16 ***
m 30.413 2.896 10.504 < 2e-16 ***
s -20.695 3.606 -5.739 3.11e-08 ***
w1 28.369 6.537 4.340 2.17e-05 ***
w2 27.480 6.548 4.197 3.93e-05 ***
w3 25.962 6.633 3.914 0.000121 ***
w4 0.890 8.173 0.109 0.913388
w5 43.524 13.567 3.208 0.001535 **
w6 41.223 14.429 2.857 0.004686 **
w7 6.091 7.262 0.839 0.402466
w8 26.610 13.604 1.956 0.051725 .
w9 39.495 14.009 2.819 0.005252 **
w10 2.265 7.551 0.300 0.764455
w11 66.661 11.478 5.808 2.19e-08 ***
w12 38.222 13.031 2.933 0.003708 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 15.62 on 221 degrees of freedom
Number of iterations to convergence: 6
Achieved convergence tolerance: 7.269e-06
# 4.1.1 - Base linear mixed-effect model ----
l.me <- lme(lfmc ~ time * leaf.type, random = ~1 | plot, data = lfmc)
# Residual checking:
plot(l.me) # heterocedasticity in the residuals
plot(ACF(l.me, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) # temporal correlation is observed at lag 1
# 4.1.2 - Variance modelling ----
l.me.vm <- update(l.me, weights = varIdent(form=~ 1 | leaf.type))
# Residual checking:
plot(l.me.vm) # variances there seem to be well modeled
AIC(l.me, l.me.vm) # modeling the variance improves the model fit
# 4.1.3 - Temporal correlation modelling ----
l.me.vm.tm <- update(l.me.vm, correlation=corARMA(p=1,q=1, form=~1))
plot(ACF(l.me.vm.tm, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) # the temporal dependence is modeled
AIC(l.me.vm, l.me.vm.tm) # modeling the temporal correlation improves the model fit
plot(l.me.vm.tm)
plot(ACF(l.me.vm.tm, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE)
hist(resid(l.me.vm.tm, type="normalized"), main="", xlab="Residuals")
qqnorm(resid(l.me.vm.tm, type="normalized"))
qqline(resid(l.me.vm.tm, type="normalized"))
# 4.1.4 - Evaluation of the fixed-effects ----
l.me.ml <- update(l.me.vm.tm, method ="ML") # the model is fitted by Maximum likelihood
l.me.ml.1 <- update(l.me.ml, ~.-time:leaf.type) # the model is reduced removing the interaction "time x leaf type"
IC_tab(l.me.ml, l.me.ml.1) # the AIC comparison suggests the interaction term to be important
M3 <- l.me.vm.tm
summary(M3)
Linear mixed-effects model fit by REML
Data: lfmc
Random effects:
Formula: ~1 | plot
(Intercept) Residual
StdDev: 3.750488 7.842661
Correlation Structure: ARMA(1,1)
Formula: ~1 | plot
Parameter estimate(s):
Phi1 Theta1
0.9076001 -0.7567060
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W M. spinosum S. bracteolactus Grass E
1.000000 2.963571 3.123740 1.070897
Fixed effects: lfmc ~ time * leaf.type
Correlation:
(Intr) time lf.tGE lf.M.s lf.S.b tm:.GE
time -0.557
leaf.typeGrass E -0.709 0.395
leaf.typeM. spinosum -0.407 0.227 0.654
leaf.typeS. bracteolactus -0.394 0.219 0.653 0.666
time:leaf.typeGrass E 0.374 -0.671 -0.592 -0.430 -0.429
time:leaf.typeM. spinosum 0.176 -0.316 -0.299 -0.742 -0.409 0.546
time:leaf.typeS. bracteolactus 0.166 -0.299 -0.318 -0.441 -0.743 0.562
t:.M.s
time
leaf.typeGrass E
leaf.typeM. spinosum
leaf.typeS. bracteolactus
time:leaf.typeGrass E
time:leaf.typeM. spinosum
time:leaf.typeS. bracteolactus 0.567
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.3449313 -0.6438200 -0.1011963 0.6158869 3.3469739
Number of Observations: 247
Number of Groups: 6
M4 <- lm(lfmc ~ time * group, data = lfmc)
summary(M4)
Call:
lm(formula = lfmc ~ time * group, data = lfmc)
Residuals:
Min 1Q Median 3Q Max
-49.673 -7.042 0.466 7.002 66.332
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 48.79275 6.44309 7.573 9.60e-13 ***
time -0.24478 0.13305 -1.840 0.06714 .
groupGW plot 5 0.49817 9.11190 0.055 0.95645
groupGW plot 6 6.14984 9.11190 0.675 0.50042
groupGE plot 1 40.32675 9.49655 4.246 3.19e-05 ***
groupSM plot 1 212.68278 9.11190 23.341 < 2e-16 ***
groupSS plot 1 227.42468 9.11190 24.959 < 2e-16 ***
groupGE plot 2 14.39850 9.49655 1.516 0.13089
groupSM plot 2 195.54780 9.11190 21.461 < 2e-16 ***
groupSS plot 2 217.12857 9.11190 23.829 < 2e-16 ***
groupGE plot 3 21.68513 9.49655 2.283 0.02334 *
groupSM plot 3 189.79613 9.49655 19.986 < 2e-16 ***
groupSS plot 3 192.53490 9.49655 20.274 < 2e-16 ***
time:groupGW plot 5 -0.01905 0.18816 -0.101 0.91944
time:groupGW plot 6 -0.10231 0.18816 -0.544 0.58718
time:groupGE plot 1 -0.80129 0.19357 -4.139 4.94e-05 ***
time:groupSM plot 1 -2.36256 0.18816 -12.556 < 2e-16 ***
time:groupSS plot 1 -2.55767 0.18816 -13.593 < 2e-16 ***
time:groupGE plot 2 -0.43459 0.19357 -2.245 0.02574 *
time:groupSM plot 2 -2.34722 0.18816 -12.474 < 2e-16 ***
time:groupSS plot 2 -2.45552 0.18816 -13.050 < 2e-16 ***
time:groupGE plot 3 -0.56657 0.19357 -2.927 0.00378 **
time:groupSM plot 3 -1.82099 0.19357 -9.407 < 2e-16 ***
time:groupSS plot 3 -2.16847 0.19357 -11.202 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 16.51 on 223 degrees of freedom
Multiple R-squared: 0.9586, Adjusted R-squared: 0.9544
F-statistic: 224.7 on 23 and 223 DF, p-value: < 2.2e-16
# Residual checking:
par(mfrow=c(2,2))
plot(M4) # a clear pattern in the residuals is observed (model assumptions are not met)
# 6 - NULL MODEL (M5) ######################################################################################
M5 <- lm(lfmc ~ 1, data = lfmc)
summary(M5)
Call:
lm(formula = lfmc ~ 1, data = lfmc)
Residuals:
Min 1Q Median 3Q Max
-88.75 -58.05 -31.45 52.66 231.06
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 95.645 4.919 19.45 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 77.3 on 246 degrees of freedom
IC_tab(M2, M3, M4, M5)#, weights = T, delta = TRUE, base=T, sort = TRUE) # the nonlinear mixed-effects model is clearly the best
IC_tab(M1)
par(mfrow=c(2,2), cex = 0.65) # (Figure 5)
v1 <- c(1, 2, 3, 4)
v2 <- c("GW", "GE", "SM", "SS")
# Nonlinear mixed-effects model (M1)
rM1 <- resid(M1, type = "normalized")
fM1 <- fitted(M1)
plot(fM1, rM1, xlab="Fitted LFMC (%)", ylab="Residuals", ylim=c(-3,3))
abline(a=0, b=0, lwd=2)
boxplot(rM1 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n", ylim=c(-3,3))
axis(side=1, at=v1, labels=v2, las=1)
plot(rM1 ~ lfmc$time, ylab="Residuals", xlab="Time (days)", ylim=c(-3,3))
abline(a=0,b=0, lw=2)
qqnorm(rM1, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM1)
# Nonlinear fixed-effects model (M2)
par(mfrow=c(2,2), cex = 0.65) # (Figure 5)
rM2 <- resid(M2)
fM2 <- fitted(M2)
plot(fM2, rM2, xlab="Fitted LFMC (%)", ylab="Residuals")
abline(a=0, b=0, lwd=2)
boxplot(rM2 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n")
axis(side=1, at=v1, labels=v2, las=1)
plot(rM2 ~ lfmc$time, ylab="Residuals", xlab="Time (days)")
abline(a=0,b=0, lw=2)
qqnorm(rM2, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM2)
par(mfrow=c(2,2), cex = 0.65)
# Linear mixed-effects model (M3)
rM3 <- resid(M3, type = "normalized")
fM3 <- fitted(M3)
plot(fM3, rM3, xlab="Fitted LFMC (%)", ylab="Residuals", ylim=c(-3,3))
abline(a=0, b=0, lwd=2)
boxplot(rM3 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n", ylim=c(-3,3))
axis(side=1, at=v1, labels=v2, las=1)
plot(rM3 ~ lfmc$time, ylab="Residuals", xlab="Time (days)", ylim=c(-3,3))
abline(a=0,b=0, lw=2)
qqnorm(rM3, main="", ylab="Sample quantiles", xlab="Theoretical quantiles")
qqline(rM3)
par(mfrow=c(2,2), cex = 0.65)
# Clasical regression (M4)
rM4 <- resid(M4)
fM4 <- fitted(M4)
plot(fM4, rM4, xlab="Fitted LFMC (%)", ylab="Residuals")
abline(a=0, b=0, lwd=2)
boxplot(rM4 ~ lfmc$leaf.type, ylab="Residuals", xlab="", xaxt="n")
axis(side=1, at=v1, labels=v2, las=1)
plot(rM4 ~ lfmc$time, ylab="Residuals", xlab="Time (days)")
abline(a=0,b=0, lw=2)
qqnorm(rM4, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM4)